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A  computational  approach  to  homogeneous  nucleation  is  proposed  based  on  Eulerian  description  of 
the  gas  phase  expansion  coupled  with  a  Lagrangian  approach  to  the  cluster  formation.  A  continuum, 
Euler/Navier-Stokes  solver  versatile  advection  code  is  used  to  model  the  gas  transport,  and  a  kinetic 
particle  solver  is  developed  in  this  work  to  simulate  cluster  nucleation  and  growth.  Parameters  in 
the  new  model  were  adjusted  so  as  to  match  the  known  theoretical  dimer  formation  equilibrium 
constants  for  the  two  gases  under  consideration,  argon  and  water.  Reasonable  agreement  between 
computed  and  available  experimental  data  was  found  in  terminal  cluster  size  distributions  for  nozzle 
water  expansions  in  a  wide  range  of  stagnation  pressures.  The  proposed  approach  was  found  to  be 
orders  of  magnitude  faster  than  a  comparable  approach  based  on  the  direct  simulation  Monte  Carlo 
method.  ©  2011  American  Institute  of  Physics,  [doi:  10.1063/1.3562370] 


I.  INTRODUCTION 

Homogeneous  condensation  plays  an  important  role 
in  many  atmospheric  and  technological  processes,  and 
understanding  its  physical  mechanisms  and  dependencies  is 
critical  for  a  number  of  engineering  applications.  One  such 
application,  pertaining  to  post  boost  vehicle  operations  at 
very  high  altitudes,  is  related  to  thruster  plume  expansion 
into  the  surrounding  rarefied  atmosphere.1  It  is  well  known 
that  particulates  of  different  kinds  are  the  main  contributor  to 
the  scattering  of  sunlight  observed  in  high  altitude  plumes. 
The  effect  of  sunlight  scattering  in  plumes  in  which  neither 
carbon  soot  nor  alumina  particles  were  present  in  significant 
concentrations,  with  the  specific  example  of  the  Apollo  8 
translunar  injection  bum,2  indicates  that  particles  must  be 
formed  in  the  rapid  expansion  of  the  exhaust  to  the  rarefied 
atmosphere,  mostly  from  the  condensation  of  water  vapor 
and  other  combustion  products  in  the  plume. 

Condensation  in  rapidly  expanding  flows  has  been  ob¬ 
served  experimentally  as  early  as  the  mid- 1930s  (Ref.  3)  and 
has  been  extensively  studied  in  the  following  decades  (see  for 
example  Ref.  4  and  the  references  therein).  Computational 
modeling  of  expanding  condensing  flows  has  a  shorter, 
although  still  a  respectable  history.  In  the  past,  two  different 
approaches  have  been  used  to  describe  homogeneous  con¬ 
densation  and,  in  particular,  cluster  nucleation  (formation 
of  small  clusters  from  monomers)  in  the  nonequilibrium 
environment  of  rapid  expansions.  In  the  first  approach,  based 
on  the  classical  nucleation  theory  (CNT)  (Ref.  5)  and  equi¬ 
librium  thermodynamics,  the  key  process  is  the  formation 
of  the  smallest  stable  droplets  possible,  so-called  critical 
clusters,  through  unimolecular  reactions  of  cluster  growth 
and  decay.  The  classical  theory  calculates  the  condensation 
and  evaporation  rates  using  the  Gibbs  distributions  and  the 
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principle  of  detailed  balance,  and  the  nucleation  rate  is  then 
calculated  assuming  a  steady  state  condition.6 

The  main  principles  of  the  classical  nucleation  theory 
in  combination  with  the  conventional  compressible  Euler 
or  Navier-Stokes  gas  dynamic  equations  were  used  by  a 
number  of  researchers  to  predict  numerically  multidimen¬ 
sional  condensing  flows  (see,  for  example,  Refs.  7-9).  The 
important  part  of  these  models  is  the  creation  of  cluster  nuclei 
at  some  critical  size  that  depends  on  local  gas  conditions. 
The  nucleation  rate  is  governed  by  CNT  and  droplet  growth 
can  be  derived  on  the  basis  of  heat  transfer  conditions 
surrounding  the  droplet  (the  description  of  Ref.  10  was  used 
in  Ref.  8).  Both  Eulerian7  and  Lagrangian8  descriptions  of 
condensed  droplets  have  been  used  in  the  literature. 

An  alternative  approach  to  modeling  homogeneous  con¬ 
densation  is  based  on  some  assumed  shape  of  the  droplet  size 
distribution,  usually  lognormal.  In  Ref.  11,  this  assumption 
is  coupled  with  a  modified  form  of  the  Hertz-Knudsen  equa¬ 
tion,  which  gives  the  droplet-gas  mass  transfer  rate  as  the 
difference  between  incoming  fluxes  from  the  gas  phase  and 
evaporative  fluxes  from  the  droplet;  a  standard  Eulerian  de¬ 
scription  was  used  to  model  the  two-phase  flow.  In  Ref.  12, 
viscous  compressible  reduced  Navier-Stokes  equations13  are 
used  for  the  gas  phase,  while  the  polydisperse  particle  behav¬ 
ior  is  described  by  the  Eulerian  aerosol  moment  model  which 
accounts  for  particle  transport  due  to  convection,  diffusion, 
inertia,  and  thermophoresis,  as  well  as  particle  dynamics  due 
to  coagulation,  nucleation,  and  condensation.  Yet  another  nu¬ 
merical  approach,  which  uses  many  of  the  CNT  assumptions 
and  has  been  applied  mostly  to  turbulent  condensing  flows,  is 
based  on  a  semi-Lagrangian  treatment  of  droplets.14,15  Semi- 
Lagrangian  methods  combine  both  Eulerian  and  Lagrangian 
points  of  view:  a  scalar  held  is  discretized  on  the  Eulerian  grid 
but  is  advanced  in  time  using  a  Lagrangian  technique. 

While  different  methods  were  applied  to  predict  cluster 
nucleation  and  growth  in  gas  hows,  most  of  the  researchers 
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that  used  the  classical  nucleation  theory  applied  the  Eulerian 
approach  to  the  gaseous  phase,  usually  based  on  the  solu¬ 
tion  of  full  or  reduced  Navier-Stokes  equations.  A  different 
strategy  was  proposed  in  Ref.  16,  where  a  particle-based  di¬ 
rect  simulation  Monte  Carlo  (DSMC)  method17  was  used  to 
compute  the  gas  flow.  A  Lagrangian  technique  was  applied  to 
model  cluster  evolution.  Similar  to  Ref.  7,  new  clusters  were 
created  at  a  critical  size,  and  their  further  growth  was  calcu¬ 
lated  with  the  CNT  approach. 

Some  important  assumptions  of  the  CNT,  such  as  uni- 
molecular  reactions  of  cluster  growth  and  decay,  and  the  use 
of  the  principle  of  detailed  balance  that  implies  thermody¬ 
namic  equilibrium,  limit  its  applicability  as  a  prediction  tool 
for  highly  nonequilibrium  flows,  such  as  rapidly  expanding 
plumes.  In  such  flows,  the  impact  of  thermal  nonequilibrium 
between  gas  and  particles  is  expected  to  significantly  impact 
the  growth  rates  and  cluster  size  distributions.  Moreover, 
the  cluster  size  distribution  may  have  a  significantly  more 
complex  shape  than  the  lognormal  distribution  often  used 
in  the  literature.  An  illustrative  example  of  such  complex 
distributions  was  provided  in  Ref.  18,  where  the  terminal 
cluster  size  distributions  were  measured  in  water  and  am¬ 
monia  expansions  for  a  wide  range  of  stagnation  pressures 
and  temperatures;  the  results  were  obtained  by  doping  the 
water  and  ammonia  clusters  by  one  Na  atom,  which  was 
photoionized  close  to  the  threshold  without  fragmentation. 

The  experimental  study18  showed  that,  for  lower 
pressures,  the  size  distribution  is  exponential;  for  higher  pres¬ 
sures,  the  size  distribution  approaches  the  lognormal  profile, 
and,  for  intermediate  pressures,  it  has  a  complex  bimodal 
shape.  The  transition  from  the  exponential  to  the  bimodal 
shape  was  explained  by  changing  governing  mechanisms 
of  cluster  growth.  For  lower  pressures,  the  clusters  grow 
mostly  through  monomer  sticking,  while  at  higher  pressures, 
the  main  mechanism  is  cluster-cluster  collisions  and  coa¬ 
lescence.  The  bimodal  shape  of  the  cluster  size  distribution 
function  for  intermediate  plenum  pressures  was  attributed 
in  Ref.  18  to  the  coalescence  of  small  particles  (such  as 
dimers  and  trimers)  on  larger  clusters  and  the  coagulation  of 
larger  clusters.  A  bimodal  distribution  of  cluster  sizes  was 
measured18  for  a  number  of  chamber  pressures,  varied  by  up 
to  an  order  of  magnitude;  typically,  it  was  observed  when 
the  average  cluster  size  was  from  below  100  to  about  1000. 
These  cluster  sizes  are  believed  to  be  largely  occurring  in  a 
number  of  applications,  including  rocket  thruster  plumes. 

The  inability  of  CNT-based  methods  to  accurately  predict 
the  cluster  size  distributions  in  strongly  nonequilibrium  flows 
dictates  the  use  of  the  second  approach,  known  as  the  kinetic 
approach,  which  treats  nucleation  as  the  process  of  kinetic 
chemical  aggregation.19  Unlike  CNT,  the  kinetic  approach 
does  not  assume  local  thermodynamic  equilibrium.  Instead, 
a  microscopic  process  of  the  interactions  of  monomers  and 
clusters  is  described  either  analytically  via  a  mathematical 
model,  e.g.,  by  the  Smoluchowski  equations  where  the  inter¬ 
action  between  particles  is  modeled  by  the  reaction  rates,20,21 
or  in  computer  simulations,  e.g.,  in  molecular  dynamics  (MD) 
calculations  where  the  interaction  is  modeled  by  an  interac¬ 
tion  potential.22,23  It  is  well  known  that  the  application  of  ei¬ 
ther  the  Smoluchowski  equations  or  the  molecular  dynamics 


approach  to  the  modeling  of  cluster  evolution  in  multidimen¬ 
sional  thruster  plume  flows  is  computationally  unfeasible. 

A  more  promising  direction  in  modeling  rapidly  expand¬ 
ing  condensing  flows  is  the  use  of  the  DSMC  method.  As 
a  numerical  approach  to  the  Boltzmann  equation,  it  is  ap¬ 
plicable  to  a  large  range  of  flow  conditions.  In  this  method, 
cluster-cluster  and  cluster-monomer  interactions  including 
the  multibody  reactions  of  cluster  nucleation  can  be  seam¬ 
lessly  incorporated.  Over  the  last  several  years,  the  DSMC 
method  has  been  extensively  and  successfully  applied  to  mod¬ 
eling  cluster  formation  and  evolution  in  supersonic  jets.24,25 
The  work  of  these  authors26  extended  the  kinetic  dimer  for¬ 
mation  approach  of  Ref.  27,  which  assumed  that  a  ternary 
collision  always  results  in  a  dimer  formation,  to  include  MD 
simulations  for  obtaining  information  on  the  probability  of 
dimer  formation  in  such  ternary  collisions.  The  work28  used 
a  temperature-dependent  probability  of  formation  of  argon 
dimers.  Another  DSMC-based  model,  which  treats  both  clus¬ 
ter  nucleation  and  evaporation  [Rice,  Ramsperger,  and  Kassel, 
or  RRK  (Ref.  29)  technique  was  used  for  the  latter]  from  the 
principles  of  the  kinetic  theory,  was  introduced  in  Ref.  30. 

The  difficulty  of  using  the  DSMC  method  as  the  mod¬ 
eling  approach  for  condensing  plumes  is  its  high  computa¬ 
tional  cost.  It  may  be  applied  to  relatively  low  density  plumes, 
when  the  typical  size  of  clusters  does  not  exceed  100-mers. 
For  higher  pressures,  this  approach  becomes  prohibitively  ex¬ 
pensive.  The  most  serious  numerical  limitation  of  the  DSMC 
method  is  related  to  the  fact  that  a  large  number  of  simulated 
particles  has  to  be  computed.  The  required  number  of  sim¬ 
ulated  particles  generally  increases  as  n2  for  2D  problems, 
and  n3  for  3D  problems,  where  n  is  the  gas  density.  Most  of 
the  simulated  particles  are  monomers;  the  statistical  scatter 
for  cluster  species  is,  therefore,  extremely  high  as  compared 
to  the  monomers.  The  use  of  species  weights  in  the  DSMC 
method  is  questionable,  since  clusters,  especially  for  larger 
pressures,  are  not  a  trace  species,  and  thus  strongly  impact  the 
flow  properties  through  the  heat  release  during  the  nucleation 
and  cluster  growth  process. 

The  main  objective  of  this  work  is  the  development  of  a 
new  method  that  would  combine  the  computational  efficiency 
of  the  Eulerian  continuum  approach  and  the  physical  accuracy 
of  the  Eagrangian  kinetic  approach.  The  proposed  method 
integrates  the  Eulerian  approach  for  monomer  gas  flow  based 
on  the  solution  of  Euler/Navier-Stokes  equations,  with  the 
Lagrangian  approach  for  cluster  formation  and  evolution 
based  on  a  DSMC-like  particle-based  algorithm.  The  work 
is  built  on  the  previous  effort30  where  the  first-principles 
model  of  homogeneous  condensation  was  formulated,  and 
all  of  the  most  important  processes  of  cluster  nucleation 
and  evolution  were  considered  at  the  microscopic  level. 
The  processes  included  in  the  model30  are  (i)  creation 
of  dimers  through  the  collision  stabilization  of  collision 
complexes,  (ii)  elastic  monomer-cluster  collisions  that 
change  the  translational  and  internal  energies  of  colliding 
particles,  (iii)  inelastic  monomer-cluster  collisions  that 
result  in  monomer  sticking,  (iv)  cluster-cluster  coales¬ 
cence,  and  (v)  evaporation  of  monomers  from  clusters. 
All  these  processes  are  present  in  the  new  method.  In 
Sec.  II,  the  details  of  the  method  are  discussed,  and  the 
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homogeneous  nucleation  rates  in  argon  and  water  thermal 
bath  environments  are  analyzed,  followed  by  the  validation 
study  that  focuses  on  comparison  of  cluster  growth  in  plumes 
with  available  experimental  data  on  terminal  cluster  size 
distribution. 

II.  NUMERICAL  APPROACH 

The  main  idea  of  the  present  numerical  method  is  to  cal¬ 
culate  gas  flow  solving  the  compressible  Euler  or  Navier- 
Stokes  equation,  model  the  nucleation  process  starting  from 
the  dimer  formation  and  up  using  the  elementary  kinetic  the¬ 
ory  for  cluster-cluster  and  cluster-monomer  collisions,  and 
exchange  the  information  between  the  continuum  and  kinetic 
parts  of  the  simulation  through  source  terms,  so  that  these 
parts  are  fully  coupled.  Similar  to  the  DSMC  method,  a  fi¬ 
nite  number  of  simulated  clusters  replace  the  real  ones,  so 
that  each  simulated  cluster  represents  a  large  number  of  real 
particles.  The  obvious  benefit  of  this  approach  is  related  to 
the  combination  of  the  numerical  efficiency  of  computing  the 
carrier  gas  flow  with  a  continuum  Eulerian  method  and  phys¬ 
ical  accuracy  of  modeling  the  cluster  evolution  with  a  kinetic 
theory  based  technique. 

A.  Eulerian  approach  to  gas  phase  expansion 

The  Eulerian-Lagrangian  approach  with  a  two-way  cou¬ 
pling  developed  in  Ref.  31  to  model  two-phase  plume  flows 
represents  the  computational  framework  of  the  new  conden¬ 
sation  model.  The  gas  flow  is  computed  using  the  Eulerian 
approach  based  on  the  solution  of  the  Euler  equations  with 
appropriate  source  terms  that  take  into  account  the  impact  of 
condensation  process  and  clusters  on  the  gas  flow: 

dp 

-£  +  v-py  =  M,  (i) 

ot 

d  px 

—  +  v  •  PW  +  X7P  =  D,  (2) 

ot 

^7  +  V  •  (e  +  p)y  =  Q,  (3) 

where  M,  D,  and  Q  are  the  corresponding  mass,  momentum, 
and  energy  source  terms  that  define  the  impact  of  condensa¬ 
tion  on  gas  molecules,  and  p,  v,  p,  and  e  are  the  gas  density, 
velocity,  pressure,  and  energy,  respectively. 

The  Euler  equations  are  solved  using  the  versatile  ad- 
vection  code  (VAC)  (Ref.  32)  modified  to  include  the 
above  source  terms.  Particle  properties  are  determined  by 
Lagrangian  tracking  of  particles  through  the  gas  flowfield 
and  statistical  averaging  of  particle  parameters.  For  the 
gas  phase,  an  explicit  time  integration  is  used,  and  the 
gas  dynamic  equations  are  solved  using  the  total  varia¬ 
tion  diminishing  (TVD)-Lax-Friedrichs  scheme  with  min- 
mod  limiter.  For  the  particle  phase,  a  fourth  order  Adams- 
Moulton  method  is  used  to  integrate  particle  equations  of 
motion. 

In  the  current  implementation,  the  clusters  are  assumed 
to  be  in  translational  equilibrium  with  the  gas,  that  is,  their 
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macroscopic  velocity  and  translational  temperature  are  as¬ 
sumed  to  be  equal  to  the  corresponding  parameters  of  the 
gas.  This  significantly  simplifies  the  implementation,  since 
then  there  is  no  need  for  the  calculation  of  cluster  drag.  At 
the  same  time  it  results  in  an  approximate  treatment  of  the 
gas-cluster  heat  and  momentum  transfer,  since  the  drag  com¬ 
ponent  is  not  included  in  the  heat  transfer  process.  Note  that 
since  the  clusters  are  small,  the  drag  contribution  to  the  heat 
transfer  is  expected  to  be  noticeably  smaller  than  that  com¬ 
ing  from  the  temperature  difference.  Note  also  that  a  more 
accurate  model  could  be  developed  that  would  take  into  ac¬ 
count  the  separation  between  cluster  and  gas  velocities  and 
the  drag  force  from  gas  to  droplets,  as  well  as  viscous  terms 
in  the  gas  phase  equations.  However,  it  is  out  of  the  scope  of 
this  work,  where  the  primary  scope  is  the  introduction  of  a 
new  particle-based  condensation  model  capable  of  a  two-way 
coupling  with  a  continuum  solver. 

At  each  time  step  At ,  the  clusters  are  moved  by  XiAt, 
where  v*  is  the  velocity  of  the  i th  cluster.  Then,  the  cluster 
collision  relaxation  processes  are  modeled  at  the  kinetic  level. 
These  processes,  that  include  the  formation  of  new  dimers, 
monomer-cluster  collisions  that  involve  energy  transfer  be¬ 
tween  internal  and  translational  modes  of  colliders,  cluster- 
cluster  coalescence,  and  cluster  growth  and  shrinking  due 
to  monomer  sticking  and  evaporation,  are  described  in  de¬ 
tail  below.  After  the  cluster  relaxation  processes,  the  changes 
in  cluster  mass  and  internal  energy  are  evaluated,  and  then 
used  to  calculate  the  right-hand  side  of  the  Navier-Stokes 
equations. 

B.  Lagrangian  approach  to  cluster  formation  and 
evolution 

Replacing  the  kinetic  modeling  of  gas  transport  with  a 
continuum  approach  is  justified  by  the  proximity  to  equilib¬ 
rium  of  velocity  distribution  functions  of  gas  molecules  in 
condensing  plumes,  where  the  gas  density  is  fairly  large,  and 
the  mean  free  path  is  many  orders  of  magnitude  smaller  than 
the  characteristic  flow  size  (usually  nozzle  throat  or  exit  di¬ 
ameter).  The  cluster  nucleation  and  evaporation  processes, 
though,  require  kinetic  treatment  for  a  number  of  reasons, 
most  notably  nonequilibrium  cluster  size  distribution  and  the 
departure  from  equilibrium  of  cluster  internal  energies.  Such 
a  kinetic,  Lagrangian  treatment  is,  therefore,  proposed  in  the 
present  work.  Although  the  kinetic,  Lagrangian  approach  to 
modeling  cluster  nucleation  follows  to  some  extent  the  first 
principle,  fully  kinetic  approach  of  the  previous  work,30  it 
has  a  number  of  key  differences,  mostly  related  to  the  fact 
that  monomers  are  simulated  at  the  continuum  level.  Thus, 
some  approximation  has  to  be  used  to  include  cluster-related 
collisions  that  involve  monomers. 


1.  Dimer  formation 

Dimers  are  formed  as  a  result  of  a  collisional  stabi¬ 
lization  of  collision  complexes  consisting  of  two  monomers 
that  collide  with  third  particles  during  their  lifetime;  the 
third  particle  is  needed  to  carry  away  extra  energy  and, 
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thus,  stabilize  the  dimer.  In  each  stabilization  event,  there  is 
also  an  energy  release  from  the  intermolecular  potential  of 
the  two  monomers  that  formed  the  collision  complex  into 
the  internal  energy  modes  of  the  newly  created  dimer  and 
the  third  particle,  and  the  translational  modes  of  their  relative 
motion. 

The  change  in  the  dimer  number  density  A  rid  as  a  result 
of  new  dimers  formed  in  each  cell  over  a  single  time  step  At 
is  calculated  from  the  known  recombination  rate  and  macro¬ 
scopic  gas  properties  in  the  cell  as  follows: 


Arid  =  KrGCn3  At. 


(4) 


Here,  KYQC  is  the  recombination  rate  constant  and  n  is  the  gas 

number  density.  Then,  the  number  of  newly  created  dimers 
is  given  by  AridVc/  Fnum,  where  Fnum  is  the  number  of  real 
clusters  represented  by  one  simulated  cluster  (similar  to  Fnum 
traditionally  used  in  the  DSMC  method)  and  Vc  is  the  cell  vol¬ 
ume.  Generally,  any  form  of  temperature  dependence  may  be 
used  to  define  ATrec;  in  this  work,  a  temperature  dependence 
similar  to  the  well  known  Arrhenius  dependence  is  used, 
Kvgc(T)  =  A  x  Bt  exp(— CT).  In  this  equation,  constants  A, 
B ,  and  C  may  be  chosen  either  from  values  known  in  the 
literature  or  selected  to  reproduce  analytical  dimer  formation 
equilibrium  constants.  Note  that  the  reverse  process  of  dimer 
loss  is  not  shown  in  Eq.  (4)  since  it  is  included  separately 
through  collisional  and  evaporative  mechanisms  discussed  in 
Secs.  II  B  3-II  B  4. 

The  initial  position  of  each  formed  dimer  is  selected  ran¬ 
domly  within  the  cell,  its  initial  velocity  is  set  equal  to  the 
macroscopic  velocity  of  the  gas  in  the  cell,  and  the  initial  clus¬ 
ter  internal  energy  is  sampled  as  follows.  First,  the  total  avail¬ 
able  energy  in  the  collision  complex-third  particle  collision  is 
assumed  to  be  equal  to 


(3  4  —  2an 

~  (2linU  + 


+ 


4  2a  c 


kT , 


(5) 


where  §int  is  the  number  of  internal  degrees  of  freedom  of 
the  monomer  (zero  if  atom),  k  is  the  Boltzmann  constant,  and 
a  is  the  variable  hard  sphere  (VHS)  model17  parameter.  Al¬ 
though  a  more  sophisticated  model  of  the  intermolecular  po¬ 
tential  may  be  used  here,  the  VHS  model  is  chosen  as  a  sim¬ 
ple  potential  that  was  proved  to  be  accurate  and  efficient  on 
kinetic  modeling  of  plume  flows.  In  this  work,  standard  VHS 
values  for  argon17  were  used,  am—ac  —  0.31.  For  water,  the 
values  am  =  ac  =  0.5  were  chosen  that  provide  a  reasonable 
approximation  of  the  viscosity-temperature  dependence  in  the 
range  of  temperatures  of  interest. 

In  Eq.  (5),  the  4  —  2am  term  corresponds  to  the  number 
of  relative  translational  degrees  of  freedom  in  the  monomer- 
monomer  collision,  and  4  —  2ac  is  the  number  of  relative 
translational  degrees  of  freedom  in  the  interaction  of  a  col¬ 
lision  complex  and  a  monomer.  After  that,  the  total  en¬ 
ergy  is  increased  by  evaporation  (dimer  dissociation)  en¬ 
ergy  F’evap  and  then  split  between  the  newly  created  dimer 
and  the  third  particle  using  the  Earsen-Borgnakke  (EB) 
(Ref.  33)  procedure.  This  procedure,  initially  developed  to 
model  energy  transfer  between  translational  and  rotational 
modes  of  colliding  molecules,  is  based  on  the  assumption  that 


after-collision  relative  translational  and  internal  energy  modes 
will  be  populated  according  to  the  local  equilibrium  distribu¬ 
tion  functions. 

In  this  work,  the  energy  transfer  was  assumed  to  include 
all  available  after-collision  energy  modes,  i.e.,  energy  of  rel¬ 
ative  motion  of  the  dimer-monomer  pair,  the  internal  energy 
of  the  monomer,  and  the  internal  energy  of  the  dimer.  The 
energy  Etot  +  £evap  is,  thus,  split  between  those  modes  using 
the  LB  approach.  For  monomers,  only  rotational  modes  are 
assumed  to  be  excited,  since  at  low  gas  temperatures  in  ex¬ 
panding  plumes  (below  300  K  in  this  work)  the  excitation  of 
vibrational  modes  is  negligible.  The  number  of  internal  de¬ 
grees  of  freedom  of  the  dimer  is  calculated  from  the  dimer 
heat  capacity  Cv  using  a  general  expression 

ru  =  /^-3,  (6) 

k 

where  i  is  the  number  of  monomers  in  the  clusters  (i  =  2  for 
dimer).  Note  that  this  expression  is  also  utilized  for  larger 
clusters.  The  values  of  heat  capacities  used  in  this  work  for 
argon  and  water  clusters  are  taken  from  Ref.  30  and  are  highly 
uncertain. 

The  dimer  formation  procedure,  thus,  results  in  the  for¬ 
mation  of  dimers  at  a  given  temperature-dependent  rate,  and 
each  of  these  dimers  is  characterized  by  a  unique  internal 
energy  that  is  subsequently  used  in  the  cluster  collision  and 
evaporation  processes. 

2 .  Inelastic  monomer-cluster  collisions 

The  interaction  between  monomers  and  clusters  is  an  im¬ 
portant  process  that  results  in  energy  transfer  between  the  in¬ 
ternal  modes  of  clusters  and  translational  modes  of  collid¬ 
ers.  The  change  in  the  cluster  internal  energies  has  a  great 
impact  on  the  evaporation  rates  and,  thus,  needs  to  be  mod¬ 
eled  properly.  It  was  pointed  out  in  Ref.  30  that  when  the  LB 
model  is  used  to  simulate  the  energy  transfer  in  monomer- 
cluster  collisions,  it  is  reasonable  to  introduce  an  inelastic 
collision  relaxation  number  Z,  which  defines  the  probabil¬ 
ity  Rinei  that  a  cluster  will  experience  an  inelastic  collision 
leading  to  a  change  in  its  internal  energy  in  a  single  collision 
as  Rinei  =  1/Z.  This  means  that  in  the  model,  only  one  out 
of  every  Z  collisions  of  a  cluster  will  change  its  internal  en¬ 
ergy,  and  all  other  collisions  will  be  elastic.  In  every  collision 
that  involves  such  a  change,  the  after-collision  energies  are 
selected  according  to  the  local  equilibrium  distribution  func¬ 
tions.  This  is  a  simplification  of  the  actual  process,  where  in¬ 
ternal  energy  change  may  occur  in  every  collision.  It  is  also 
similar  to  the  rotational  and  vibrational  relaxation  numbers 
Zr  and  Zv  widely  used  in  the  DSMC  method.  Values  of  Z 
were  proposed  in  Ref.  30  that  provide  good  agreement  with 
known  theoretical  dimer  formation  equilibrium  constants  for 
argon  and  water;  these  values  are  used  in  this  work.  For  wa¬ 
ter,  Z  =  10  over  all  temperatures,  and  for  argon,  it  depends 
on  gas  temperature  T ,  increasing  from  3  at  T  ~  0  to  25  at  T 
=  500  K. 

A  similar  approach  is  used  in  the  present  model,  with  one 
significant  exception.  Among  many  monomer-cluster  colli¬ 
sions,  only  those  that  cause  the  cluster  internal  energy  change 
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FIG.  1.  Flow  chart  for  cluster-monomer  inelastic  collision  modeling. 


are  important  in  terms  of  cluster  evolution,  and,  thus,  only 
those  collisions  were  modeled.  Since  the  value  of  Z  is  typi¬ 
cally  noticeably  larger  than  one,  such  an  approach  allows  for 
significant  reduction  in  computational  time.  The  algorithm, 
therefore,  is  reduced  to  the  analysis  of  each  cluster  in  the 
computational  domain  in  terms  of  possible  inelastic  collisions 
with  monomers  as  follows. 

First,  note  that  the  probability  that  a  cluster  will  experi¬ 
ence  an  inelastic  collision  leading  to  a  change  in  its  internal 
energy  during  time  t  is  equal  to  p  =  1  —  exp(— vz / Z),  where 
v  is  the  collision  frequency  of  clusters  on  monomers.  There¬ 
fore,  the  time  to  the  next  inelastic  collision  can  be  sampled 
as  Tnext  =  —  Z/v  log(9^),  where  9^  is  the  random  number  uni¬ 
formly  distributed  between  0  and  1 .  The  algorithm  to  model 
inelastic  collisions  for  a  given  cluster  over  time  step  At  is 
given  in  Fig.  1. 

Unlike  Ref.  30,  where  the  kinetic,  microscopic  informa¬ 
tion  on  monomers  that  include  their  individual  energy  states 
and  velocities  is  available,  the  present  method  provides  only 
the  macroscopic  information  such  as  temperature  and  number 
density.  While  sufficient  to  calculate  the  local  collision  fre¬ 
quency  and  temperature-dependent  internal  energy  relaxation 
number  Z,  this  is  not  enough  to  simulate  monomer-cluster 
collisions  at  the  kinetic  level.  In  order  to  do  that,  additional  in¬ 
formation  about  the  velocity  and  energy  distribution  functions 
of  monomers  is  necessary.  In  this  work,  the  internal  energy  of 
the  colliding  monomer  and  the  relative  translational  energy 
of  the  colliding  monomer-cluster  pair  are  sampled  from  the 
corresponding  equilibrium  distributions.  Note  that  it  is  also 
possible  to  use  the  Chapman-Enskog  distribution  function34 
for  the  translational  modes,  but  the  results  in  most  situations 
are  not  expected  to  be  significantly  different.  The  total  col- 
lisional  energy,  which  is  the  sum  of  these  two  energies  and 
internal  energy  of  the  cluster,  is  then  redistributed  between 
the  relative  translational  and  the  internal  modes  of  the  cluster 
and  the  monomer  using  the  LB  model.  The  numbers  of  the 
corresponding  degrees  of  freedom  are  defined  as  described  in 
the  dimer  formation  section. 
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3.  Cluster  growth  and  evaporation 

The  key  processes  that  determine  small  cluster  evolution 
are  sticking  and  evaporation  of  monomers  off  the  clusters.  In 
order  to  drastically  reduce  the  requirements  to  the  minimum 
time  step  used  in  the  simulation  and  provide  an  accurate  ac¬ 
count  of  the  evaporation,  and  sticking  events  of  a  single  clus¬ 
ter,  the  growth  and  evaporation  processes  are  combined  in  a 
single  procedure  as  follows. 

The  cluster  sticking  rate  is  calculated  as  vs  =  nPs(crcg), 
where  Ps  is  the  probability  that  a  monomer  will  stick  to  the 
cluster  after  the  collision,  oc  is  the  monomer-cluster  colli¬ 
sion  cross  section  (in  this  work  it  is  calculated  using  the  hard 
sphere  model)  and  g  is  the  relative  collision  velocity.  In  the 
hard  sphere  model,  where  the  collision  cross  section  is  writ¬ 
ten  as  7T  J2,  the  collision  diameter  d  is  given  as  the  average  of 
the  diameter  of  the  colliding  monomer  obtained  from  the  VHS 
model17  (4.17  A  for  argon  and  6.2  A  for  water)  and  the  clus¬ 
ter  diameter  obtained  through  an  empirical  correlation  used 
extensively  in  the  past  (see,  for  example,  Ref.  16): 

dc  =  2x(A-i1/3  +  B),  (7) 

where  A  and  B  are  species-dependent  constants.  In  this  work, 
the  values  of  A  and  B  were  2.3  x  10“ 10  and  3.4  x  10“ 10  m 
for  argon,24  and  1.9  x  10-10  and  2.4  x  10-10  m  for  water.35 

For  the  water  monomer-cluster  sticking  collision  proba¬ 
bility,  the  empirical  dependence  of  the  probability,  e  on  the 
species  diameter  d  and  mass  m  given  in  Ref.  36  is  used.  After 
a  simple  transformation,  one  may  write 

dl  t  m  v/2  (8) 

(dt  +  di)2  \mi  +  mi  J 

where  indices  i  and  1  refer  to  the  cluster  of  size  i  and 
the  monomer,  respectively.  For  argon  monomer  sticking,  the 
size  dependent  probability  of  Ref.  30  is  used,  where  the 
sticking  probability  increases  with  cluster  size,  from  0.06 
for  monomer-dimer  collisions  to  0.9  for  monomer- 1 5 -mers 
collisions. 

To  evaluate  the  rate  of  evaporation  of  monomers  from  the 
cluster  surface,  the  RRK  model29  is  used,  similar  to  Ref.  30. 
Following  Ref.  37,  this  work  calculates  the  evaporation  rate 
ke  using 

1  at  ( ^int  _  ^evap  \ 

ke  =  vNs  - - - -  .  (9) 

\  Bint  J 

Here,  v  is  the  vibration  frequency,  Ns  is  the  number  of 
surface  atoms,  and  E-mi  is  the  cluster  internal  energy.  For 
dimers,  the  exponent  3 i  —  7  is  replaced  with  1.  The  number 
of  surface  atoms  Ns  is  i  for  i  <  5,  i  —  1  for  4  <  i  <7,  and 
(367t)1/3(/1/3  —  l)2  for  i  >  6.  The  vibration  frequency  was 
taken  to  be  1012  s_1  for  argon  clusters37  and  2.68  x  1012  s_1 
for  water  clusters.38 

With  the  evaporation  and  sticking  rates  defined  by  the 
above  expressions,  the  algorithm  used  to  model  sticking  and 
evaporation  processes  is  given  in  Fig.  2.  For  cluster  growth 
(monomer  sticking),  the  monomer  internal  energy  and  rel¬ 
ative  translational  energy  are  sampled  from  the  correspond¬ 
ing  equilibrium  distributions,  and  the  after- sticking  cluster  in¬ 
ternal  energy  is  equated  to  the  sum  of  cluster  precollisional 
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FIG.  2.  Flow  chart  for  cluster  growth  and  evaporation  modeling. 


internal  energy,  internal  energy  of  the  monomer,  relative 
translational  energy,  and  evaporation  energy  £evap- 

For  cluster  evaporation,  the  cluster  internal  energy  is  de¬ 
creased  by  £evap,  and  the  remaining  energy  is  redistributed 
between  the  cluster  internal  modes,  internal  modes  of  the  de¬ 
parting  monomer,  and  relative  translational  modes  using  the 
LB  model.  Note  here  that  only  the  cluster  internal  energy  is 
calculated,  while  the  cluster  velocity  and  monomer  properties 
are  assumed  to  accommodate  to  the  gas  properties. 

4.  Cluster-cluster  collisions 

The  cluster-cluster  collisions  must  be  taken  into  ac¬ 
count  for  an  accurate  description  of  cluster  evolution,  since 
it  is  the  key  factor  determining  the  size  distribution  of  larger 
clusters.  Cluster-cluster  collisions  have  different  outcomes, 
which  generally  may  be  classified  as  either  coalescence  or 
reflexive  and  stretching  separations.  The  dynamics  of  water 
droplet  collisions  for  macroscopic  particles  was  studied  ex¬ 
perimentally  in  Ref.  39  where  the  boundaries  between  both 
of  the  separating  collisions  and  the  coalescence  collision  were 
examined  as  a  function  of  the  size  ratio  and  the  Weber  num¬ 
ber  in  the  wide  range  of  Weber  numbers  from  1  to  100.  For 
microscopic  particles  with  sizes  from  dimers  to  1000-mers, 
the  authors  are  not  aware  of  any  comparable  to  the  Ref.  39 
systematic  study  where  the  results  of  cluster  collisions  would 
be  analyzed  for  different  Weber  numbers.  Extrapolating  the 
results  of  Ref.  39  to  microscopic  particles  of  interest  in  this 
work,  one  can  notice  that  for  typical  plume  temperatures  on 
the  order  of  100  K  and,  thus,  Weber  numbers  on  the  order  of 
the  unity  or  less,  the  clusters  would  mostly  experience  coa¬ 


lescence  and  not  separation.  However,  such  an  extrapolation, 
although  partially  justified  for  100-mers  and  1000-mers,  is 
much  more  questionable  for  smaller  clusters,  where  more  re¬ 
flexive  collisions  may  be  expected.  In  this  work,  in  the  ab¬ 
sence  of  reliable  size  and  relative  velocity  dependence  of  the 
collision  outcome  for  small  clusters,  a  constant  coalescence 
probability  is  assumed. 

The  cluster-cluster  collisions  are  modeled  using  the 
conventional  DSMC  algorithm.  The  majorant  frequency 
scheme40  of  the  DSMC  method  was  utilized  for  this  purpose. 
At  every  time  step,  the  maximum  cluster  size  is  obtained  in 
each  cell.  Then,  the  majorant  collision  frequency  is  calcu¬ 
lated  based  on  this  maximum  cluster  size  and  maximum  rel¬ 
ative  collision  velocity  evaluated  from  the  local  gas  temper¬ 
ature.  The  majorant  collision  frequency  is  then  multiplied  by 
the  coalescence  probability,  since  only  the  coalescence  events 
are  modeled  (reflexive  separation  is  believed  to  have  negligi¬ 
ble  effect  on  cluster  properties,  and  the  stretching  separation 
process  is  not  included  in  the  present  model).  After  a  pair  of 
clusters  K  and  L  is  selected  for  physical  collision,  the  coa¬ 
lescence  event  is  modeled,  with  the  result  being  a  larger  clus¬ 
ter  M  with  mass  m  and  internal  energy  E-mt  calculated  from 
the  properties  of  colliding  clusters  using  the  mass  and  energy 
conservation  constraints.  The  laws  dictated 

TYlM  —  M  K  +  Wl  L  i  =  E[nt  k  T  ^int,L  Q> 

with  Q  =  -Qm  +  Qk  +  Ql,  where  Ql  is  the  energy  of  va¬ 
porization  of  cluster  i . 

After  all  collision  and  evaporation  processes  are  simu¬ 
lated  for  a  given  time  step,  the  mass  and  energy  changes  over 
this  time  step  are  calculated  over  all  cells  in  order  to  be  in¬ 
cluded  in  the  Eulerian  gas  flow  equations.  The  primary  pur¬ 
pose  of  this  step  is  accurate  conservation  of  all  conservative 
properties  in  the  simulation. 

III.  THERMAL  BATH  RELAXATION 

Inelastic  cross  sections  for  monomer-monomer  and 
monomer-cluster  collisions  are  needed  as  part  of  a  compre¬ 
hensive  validation  of  a  kinetic  condensation  model.  These 
cross  sections,  in  general  a  function  of  the  translational  and 
internal  energy  states  of  precollision  and  postcollision  parti¬ 
cles,  are  unavailable  for  the  species  and  temperatures  desired. 
However,  equilibrium  rates  for  nucleation  and  evaporation  for 
both  water  and  argon  are  available  in  literature.  Furthermore, 
a  necessary  condition  for  the  model  to  fulfill  is  that  it  produces 
correct  behavior  in  equilibrium,  although  it  does  not  guaran¬ 
tee  the  model  will  have  correct  nonequilibrium  behavior. 

Such  equilibrium  behavior  was  modeled  in  this  work  by 
a  thermal  bath  relaxation  of  both  argon  and  water  at  vari¬ 
ous  temperature  conditions.  The  equilibrium  constants  for  the 
formation  of  clusters  were  calculated  and  compared  to  pub¬ 
lished  results  of  Refs.  41-43.  In  addition,  they  were  com¬ 
pared  to  the  previously  obtained  results  of  the  DSMC -based 
model.30 

To  model  argon  equilibrium,  over  1  x  106  simulated  par¬ 
ticles  were  used,  and  the  run  was  allowed  to  run  until  a 
steady  equilibrium  value  was  reached,  usually  about  half  a 
million  time  steps.  The  time  step  for  argon  was  selected  to  be 
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FIG.  3.  Argon  dimer  formation  equilibrium  constant  as  a  function  of  gas 
temperature. 


Temperature  (K) 


FIG.  4.  Water  dimer  formation  equilibrium  constant  as  a  function  of  gas  tem¬ 
perature. 


10“ 10  s,  so  that  there  was  on  average  fewer  than  one  collision 
per  molecule  per  time  step,  making  the  results  independent  of 
step  size.  A  monomer  number  density  of  1022  molecules/m3 
was  selected  to  ensure  that  clusters  made  up  less  than  0. 1%  of 
the  gas  while  maintaining  the  10  x  106  particle  requirement. 
This  ensured  that  adequate  numbers  of  particles  were  present 
for  statistics,  and  that  the  clusters  did  not  have  a  significant 
impact  on  the  behavior  of  the  gas. 

The  dimer  formation  rate  krec  for  argon  was  computed 
using  the  stable  dimer  formation  rate  from  Ref.  41,  which  is 
written  as 

kiec  =  Ax  Bt exp (-CT).  (10) 

The  values  of  A,  B ,  and  C  given  in  Ref.  41  are  A  =  10.15 
x  10-44m6  molecules-2  s-1,  B  =  —0.278  and  C  =  3.10 
x  10-3  K-1. 

The  argon  equilibrium  constant  as  a  function  of  temper¬ 
ature  is  shown  in  Fig.  3.  It  is  compared  with  the  DSMC  re¬ 
sults  from  Ref.  30  and  the  theoretical  predictions  of  Ref.  42. 
There  is  good  agreement  between  this  model  and  the  theo¬ 
retical  predictions,  although  it  becomes  somewhat  worse  at 
higher  temperatures.  Since  the  condensation  of  argon  is  very 
small  at  those  temperatures,  such  a  deviation  is  not  expected 
to  matter  in  plume  flows.  Note  also  that  the  present  model 
uses  the  values  of  the  inelastic  collision  number  Z  as  a  func¬ 
tion  of  temperature  taken  from  Ref.  42,  where  this  parameter 
was  adjusted  to  fit  the  equilibrium  constant  of  Ref.  42  in  the 
range  of  temperatures  between  100  and  300  K.  No  fitting  of 
any  numerical  or  physical  parameter  of  the  present  model  to 
match  the  theoretical  equilibrium  constants  was  conducted  for 
argon,  which  resulted  in  some  difference. 

In  modeling  water  equilibrium,  about  1  x  106  simulated 
particles  were  used,  and  the  process  was  again  allowed  to 
run  until  a  steady  equilibrium  was  reached,  about  half  a  mil¬ 
lion  time  steps.  The  time  step  was  set  at  1  x  10-8  s,  which 
ensured  not  only  that  there  was  still  much  fewer  than  one 
collision  per  particle  per  time  step,  but  also  that  the  system 


would  reach  complete  equilibrium  within  a  reasonable  num¬ 
ber  of  time  steps.  A  number  density  of  1021  molecules/m3 
was  selected  to  ensure  that  clusters  were  less  than  0.1%  of 
the  gas  while  maintaining  lxlO6  simulated  particles.  Due  to 
the  low  fraction  of  clusters,  this  ensured  that  the  behavior  of 
the  gas  was  not  influenced  by  the  presence  of  the  clusters. 
Note  that  using  a  smaller  gas  density  does  not  change  the  re¬ 
sults  of  the  simulations. 

The  authors  are  not  aware  of  any  theoretical,  exper¬ 
imental,  or  computational  results  that  would  provide  the 
temperature-dependent  rate  of  dimer  formation,  similar  to  that 
of  Ref.  41  for  argon.  Although  it  is  generally  possible  to  use 
a  rate  obtained  from  the  kinetic  theory  for  binary  and  ternary 
collisions,  such  a  rate  would  inherently  include  a  number  of 
assumptions  and  adjustable  parameters,  such  as  the  dimer  sta¬ 
bilization  probability,  the  collision  complex  lifetime,  and  oth¬ 
ers,  that  would  affect  the  obtained  equilibrium  constant.  It  ap¬ 
pears  more  reasonable,  therefore,  to  use  a  dimer  formation 
rate  that  provides  acceptable  agreement  with  available  dimer 
formation  equilibrium  constants.  In  this  work,  the  Arrhenius 
expression  given  by  Eq.  (10)  was  used,  with  constants  ad¬ 
justed  to  produce  the  given  equilibrium  constants:  A,  B,  and 
C  are  found  to  be  A  =  5.42  x  10-41  m6  molecules-2  s-1,  —1, 
and  2.2  x  10-3  K-i. 

The  water  equilibrium  constant  is  shown  in  Fig.  4.  It  is 
compared  with  the  theoretical  predictions43  and  the  numerical 
results  of  the  DSMC-based  model.30  As  expected  from  using 
adjusted  constants  in  Eq.  (10),  there  is  excellent  good  agree¬ 
ment  between  the  current  model  and  the  results  of  Ref.  43  at 
all  temperatures  investigated  in  this  work.  Furthermore,  the 
present  model  provides  a  much  better  ability  to  match  dimer 
formation  equilibrium  constants  than  Ref.  30. 

IV.  WATER  CLUSTER  SIZE  DISTRIBUTION  IN  NOZZLE 
EXPANSION 

The  second  part  of  the  validation  and  numerical  anal¬ 
ysis  of  the  presented  condensation  model  is  focused  on  the 
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FIG.  5.  Gas  temperature  profile  along  the  nozzle  axis  for  po  =  1.577  bars  (left)  and  po  =  8.307  bars  (right). 


nucleation  and  evolution  of  small  water  clusters  in  a  conical 
nozzle.  The  study  was  prompted  by  the  availability  of  high 
quality  experimental  data18  on  the  terminal  size  distribution 
of  water  clusters  in  the  wide  range  of  flow  conditions  where 
the  cluster  size  distribution  changes  its  shape  from  exponen¬ 
tial  at  low  pressures  to  bimodal  at  intermediate  pressure  to 
lognormal  at  high  pressures.  The  experimental  results  were 
obtained  by  doping  the  water  clusters  by  one  Na  atom,  which 
is  photoionized  close  to  the  threshold  without  fragmentation. 
The  nozzle  has  a  conical  diverging  section  with  a  41°  open¬ 
ing  angle,  a  total  length  of  2  mm,  and  a  throat  diameter  of 
50  [i m.  Four  different  stagnation  pressures  were  computed, 
considered  in  Ref.  18,  1.577,  2.173,  5.144,  and  8.307  bars, 
with  the  corresponding  stagnation  temperature  of  495  K. 
Since  the  background  pressure  effect  in  the  experiment  is  be¬ 
lieved  to  be  small,44  expansion  into  a  vacuum  is  modeled. 

The  computations  were  conducted  on  a  500  x  150  spa¬ 
tial  grid,  with  cell  sized  reduced  in  the  radial  direction  and  in¬ 
creased  in  the  axial  direction.  Previous  studies3'  have  shown 
an  insignificant  impact  of  the  wall  conditions  on  the  coreflow 
where  the  cluster  sizes  are  recorded,  so  that  the  use  of  the 
Euler  solver  for  the  gas  phase  is  reasonable.  The  number  of 
simulated  droplets  was  about  500  000,  which  was  found  to 
provide  adequate  statistical  accuracy  for  the  calculations.  The 
particles  were  assumed  to  condense  on  the  nozzle  surface. 
Uniform  inflow  conditions  were  imposed  at  the  nozzle  throat, 
calculated  from  the  isentropic  flow  relations.  To  compare 
the  cluster  size  distributions  with  the  terminal  distributions18 
measured  far  downstream  from  the  nozzle,  the  computed  size 
distributions  at  several  stations  along  the  nozzle  axis  were  an¬ 
alyzed  to  provide  truly  terminal,  distance-independent  distri¬ 
butions.  The  domain  size  was  increased  in  the  axial  direction 
from  4  mm  for  the  lowest  pressure  to  20  mm  for  the  highest 
pressure  to  ensure  that  the  size  distributions  at  the  exit  bound¬ 
ary  are  essentially  frozen. 

A  typical  run  time  for  the  lowest  pressure  under  consid¬ 
eration  was  several  hours  and  for  the  highest  pressure  was  up 
to  2  days  on  a  single  processor  computer.  Comparing  these 
numbers  with  those  of  Ref.  30  where  a  DSMC  method  was 


used  to  model  a  1.577  bars  water  expansion,  one  finds  that  the 
new  approach  is  about  50  times  faster  than  the  DSMC -based 
method  for  the  lowest  pressure.  This  factor  will  grow  signifi¬ 
cantly  with  pressure.  The  reduction  in  run  time  is  mostly  re¬ 
lated  to  the  time  efficient  modeling  of  gas  transport  with  a 
continuum  method.  Since  clusters  comprise  only  a  relatively 
small  fraction  of  the  particles  in  the  flow,  gas  transport  model¬ 
ing  is  the  most  time  consuming  part  of  any  DSMC-based  tech¬ 
nique.  Note  that  species  weights  for  cluster  species  would  re¬ 
duce  the  time  requirements  of  the  DSMC-based  condensation 
model  but  that  application  of  weights  is  questionable  in  con¬ 
densing  flows  since  the  condensation  significantly  changes 
the  gas  flow. 

Consider  first  the  gas  and  particle  properties  along  the 
nozzle  axis.  The  gas  translational  temperature  for  the  lowest 
and  highest  pressures  under  consideration  is  shown  in  Fig.  5. 
Here,  X  =  0  corresponds  to  the  nozzle  throat.  As  expected, 
the  water  nucleation  results  in  a  noticeable  increase  in  gas 
temperature  due  to  condensation.  For  the  1.577  bars  case,  the 
temperature  in  the  plume  region  is  up  to  30  K  higher  when  the 
condensation  is  included,  which  is  comparable  to  the  magni¬ 
tude  of  the  temperature  in  the  noncondensing  flow.  The  small 
temperature  increase  at  about  0.25  mm  is  related  to  the  com¬ 
pression  wave  that  originates  near  the  nozzle  throat  and  prop¬ 
agate  to  the  nozzle  axis.  It  is  present  both  in  the  condensing 
and  noncondensing  flow,  and  the  location  is  nearly  the  same 
since  the  impact  of  the  condensation  is  not  very  significant  at 
this  point. 

For  the  8.307  bars  case,  the  influence  of  the  condensa¬ 
tion  is  obvious  almost  immediately  after  the  nozzle  throat 
(the  temperatures  start  to  deviate  after  the  first  100  ji m  from 
the  throat),  and  in  the  plume  the  gas  temperature  is  several 
times  higher  in  the  condensing  flow.  Some  statistical  scatter  is 
seen  in  this  figure,  where  instantaneous  gas  properties  are  pre¬ 
sented  (an  interpolation  procedure  was  used  here  to  smooth 
the  results).  The  instantaneous  properties  are  dependent  on 
current  cluster  properties,  and  the  use  of  only  a  finite  number 
of  simulated  clusters  contributes  to  their  scatter.  The  higher 
temperatures  in  the  nozzle  for  the  condensing  flow  cause  the 
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FIG.  6.  Computed  cluster  size  distributions  for  different  coalescence  probabilities  P  at  po  =  1.577  bars  (left)  and  po  =  2.173  bars  (right). 


formation  of  a  compression  wave  near  the  nozzle  lip,  that 
propagates  downstream  and  reflects  at  the  axis  at  X  ~  12  mm. 
This  results  in  a  significant  rise  in  gas  temperature. 

Consider  now  the  terminal  cluster  size  distributions  at 
different  stagnation  pressures.  Note  first  that  there  are  several 
important  properties  that  strongly  affect  the  size  distribu¬ 
tions,  among  which  are  the  evaporation  heat,  heat  capacity, 
monomer  sticking,  and  cluster  coalescence  probabilities. 
The  first  two  of  these  properties  are  mostly  functions  of  the 
cluster  size,  and  the  latter  two,  being  characteristics  of  binary 
collisions,  depend  on  the  cluster  sizes,  internal  energies,  and 
relative  collision  velocities.  The  use  of  a  constant  coalescence 
probability  in  this  work  is  a  significant  oversimplification 
of  the  actual  cluster  collision  process,  primarily  related 
to  the  lack  of  information  on  collisions  of  small  clusters. 
While  the  coalescence  probability  of  two  relatively  large 
clusters  (100-mers  and  larger)  may  be  reasonably  assumed 
to  be  close  to  the  unity  for  Weber  numbers  on  the  order  of 
1,  the  coalescence  of  smaller  clusters  is  less  likely  and  for 


the  limiting  case  of  dimer  collisions  may  approach  that  of 
monomer  sticking,  which  is  about  0.2  for  water. 

The  numerical  analysis  has  shown  that  the  size  distribu¬ 
tion  significantly  depends  on  the  coalescence  probability,  see 
Fig.  6.  The  increase  in  the  coalescence  probability  from  0.25 
to  1  results  in  a  significant  redistribution  of  cluster  sizes  and 
a  shift  from  smaller  sizes  to  larger  ones.  Such  a  trend  is  ex¬ 
pected,  since  a  higher  coalescence  probability  at  a  given  colli¬ 
sion  rate  increases  the  population  of  large  clusters.  Although 
the  coalescence  is  accompanied  by  energy  release  from  the 
electron  structure  of  smaller  clusters  to  the  internal  energy  of 
larger  clusters,  the  larger  internal  energy  is  then  redistributed 
over  a  significantly  larger  number  of  internal  degrees  of  free¬ 
dom.  The  resulting  gas  temperatures  were,  therefore,  found 
not  to  change  noticeably  with  the  coalescence  probability. 

Comparison  of  the  computed  and  experimental  cluster 
size  distributions18  for  these  pressures  is  presented  in  Fig.  7. 
Since  the  dimers  were  not  measured  in  Ref.  18;  hereafter,  the 
experimental  points  were  normalized  to  match  the  population 
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FIG.  7.  Terminal  cluster  size  distributions  for  po  =  1.577  bars  (left)  and  po  =  2.173  bars  (right):  comparison  with  data  (Ref.  18). 
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FIG.  8.  Terminal  cluster  size  distributions  for  po  =  5.144  bars  (left)  and  po  =  8.307  bars  (right):  comparison  with  data  (Ref.  18). 


of  the  computed  clusters  excluding  dimers.  The  case  with  a 
coalescence  probability  of  0.5  gives  better  agreement  with 
the  data  and  is,  therefore,  shown  here.  It  needs  to  be  men¬ 
tioned  that  the  local  minimum  observed  for  6-mers  and  a  lo¬ 
cal  maximum  observed  for  8-mers  are  not  statistical  fluctua¬ 
tions,  but  the  consequence  of  the  corresponding  minimum  and 
maximum  in  the  cluster  evaporation  energies.  Note  that  for 
Po  =  1.577  bars,  the  best  agreement  with  the  data  would  pro¬ 
duce  a  computation  that  utilizes  a  constant  coalescence  prob¬ 
ability  between  0.25  and  0.5,  whereas  for  po  =  2.173  bars, 
a  larger  coalescence  probability  between  0.5  and  1.0  would 
produce  a  better  agreement.  This  is  reasonable,  since  higher 
pressures  are  generally  characterized  by  higher  degree  of  nu- 
cleation  and  larger  cluster  sizes,  for  which  the  coalescence 
probability  is  expected  to  increase. 

For  the  two  largest  pressures  under  consideration,  the 
computations  with  a  coalescence  probability  of  1  provide  bet¬ 
ter  agreement  with  the  data,  and  the  corresponding  results  are 
shown  in  Fig.  8.  For  po  =  5.144  bars,  the  computed  location 
of  the  second  maximum  in  the  distribution  function  agrees 
well  with  the  corresponding  experimental  value,  although  the 
population  of  such  clusters  is  somewhat  higher  in  the  experi¬ 
ments.  The  most  noticeable  difference  is  observed  in  the  large 
cluster  tail,  where  clearly  more  clusters  were  observed  in  the 
experiment.  In  the  calculation,  the  large  cluster  tail  is  closer 
to  the  lognormal  shape.  Interestingly,  the  situation  is  oppo¬ 
site  for  po  =  8.307  bars,  for  which  the  tail  is  somewhat  more 
populated  in  the  numerical  prediction.  More  importantly,  the 


TABLE  I.  Average  computed  and  measured  cluster  sizes  at  different 
pressures. 


Stagnation 
pressure  (bars) 

Computed 

Measured 

1.577 

12 

9 

2.173 

18 

20 

5.144 

107 

80 

8.307 

417 

338 

numerical  results  do  not  produce  a  clear  bimodal  structure 
at  this  pressure.  Although  this  is  clearly  related  to  some  ap¬ 
proximations  used  in  the  model,  more  research  is  needed  to 
single  out  the  most  important  reason  for  this.  The  average 
cluster  sizes  for  the  above  computation  versus  experiment 
comparisons  are  summarized  in  Table  I.  There  is  reasonable 
agreement  between  the  results,  especially  for  the  three  lowest 
pressures. 

V.  CONCLUSIONS 

A  new  method  for  modeling  homogeneous  condensation 
is  presented,  based  on  the  Eulerian  description  of  the  gas 
phase  coupled  with  the  Lagrangian  approach  to  the  cluster 
phase  formation.  A  continuum,  Euler/Navier-Stokes  solver 
VAC  is  used  to  model  the  gas  transport,  and  a  kinetic  particle 
solver  is  developed  in  this  work  to  simulate  cluster  nucleation 
and  growth.  Conservation  of  mass,  momentum,  and  energy 
is  enforced  through  a  two-way  coupling,  with  gas  properties 
influencing  the  cluster  evolution  through  the  dimer  formation 
and  monomer-cluster  collisions  (both  elastic  and  inelastic), 
and  mass,  momentum,  and  energy  transfer  from  the  cluster 
to  the  gas  phase  handled  via  source  terms  in  the  continuum 
equations.  The  proposed  approach  is  orders  of  magnitude 
faster  than  a  comparable  approach  based  on  the  DSMC 
method.  Note  also  that  it  may  easily  be  extended  to  model 
heterogeneous  condensation. 

The  following  cluster-related  processes  are  taken  into 
account  in  the  kinetic  solver:  (i)  collisional  dimer  forma¬ 
tion  that  uses  theoretical  temperature-based  dimer  formation 
rates  defining  the  number  of  dimers  created  in  each  cell  per 
time  step,  (ii)  elastic  monomer-cluster  collisions  that  change 
the  translational  and  internal  energies  of  colliding  particles, 
with  energy  transfer  modeling  using  the  Larsen-Borgnakke 
model,  (iii)  inelastic  monomer-cluster  collisions  that  result  in 
monomer  sticking,  (iv)  cluster-cluster  coalescence  simulated 
with  a  conventional  DSMC  collision  algorithm  based  on  the 
majorant  frequency  scheme,  and  (v)  evaporation  of  monomers 
from  clusters  based  on  the  RRK  model. 
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The  new  model  was  found  to  reproduce  the  well  known 
theoretical  dimer  formation  equilibrium  constants  for  the  two 
gases  under  consideration,  argon  and  water.  Water  nozzle  ex¬ 
pansion  was  modeled  with  the  stagnation  pressure  ranging 
from  1.5  to  8.3  bars,  which  corresponds  to  the  average  clus¬ 
ter  size  increasing  from  below  10  to  over  300.  The  results 
on  the  terminal  cluster  were  found  sensitive  to  the  cluster  co¬ 
alescence  probability,  with  the  average  cluster  size  increas¬ 
ing  significantly  when  this  probability  was  increased  from 
0.25  to  1.  Comparisons  with  available  experimental  data  have 
shown  good  agreement  at  lower  pressures,  and  somewhat 
worse  agreement  at  the  highest  pressure  under  consideration, 
where  no  visible  bimodal  size  distribution  structure  was  no¬ 
ticed  in  the  calculations. 
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